﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace LxCFD
{
    public partial class ZoneCFDFields
    {
        public VelocityVector Velocity;
        public PressureCorrectionBasic PCor;
        public PressureSIMPLE P;

        public void CheckCreate_Vel(bool isSupportByParent)
        {
            bool advice = isSupportByParent && PhysicalModel.PhysicalType == PhysicalTypes.Fluid;
            if (advice)
            {
                if (Velocity == null)
                {
                    Velocity = new VelocityVector();
                    Velocity.CheckCreate(Zone);
                    PCor = new PressureCorrection();
                    PCor.CheckCreate(Zone);
                    P = new PressureSIMPLE();
                    P.CheckCreate(Zone);
                }
                AppRes.SolveCenter.SIMPLEGroup.SIMPLEs.Add(PCor);
                AppRes.SolveCenter.InitializeGroup.Initialization_Normal6.Add(Velocity);
                AppRes.SolveCenter.InitializeGroup.Initialization_Normal6.Add(P);
                AppRes.MonitorCenter.Vel_U.Add(Velocity.U);
                AppRes.MonitorCenter.Vel_V.Add(Velocity.V);
                AppRes.MonitorCenter.Vel_W.Add(Velocity.W);
                AppRes.MonitorCenter.Continuity.Add(Velocity);
            }
            else
            {
                Velocity = null;
                PCor = null;
                P = null;
            }
            CreateFields_TurbBase(advice);
        }
        public void SetParameter_Vel()
        {
            if (PhysicalModel.PhysicalType != PhysicalTypes.Fluid) return;
            Value1OppositeFromArray1.Vector UVW_PSource = new Value1OppositeFromArray1.Vector();
            UVW_PSource.setParameter(P.Grad);
            var uvws = new ValueListMaker.Vector();
            uvws.AddValueProvider(UVW_PSource);
            if (PhysicalModel.IsGravityEnabled && PhysicalModel.IsEnergyEnabled && PhysicalModel.IsBoussinesqEnabled)
            {
                if (Density is Value1Const.Double vd)
                {
                    var value3D1Base = new BoussinesqSource(PhysicalModel.ReferenceTemperature, Material.ThermalExpansionCoefficent, vd.Number, PhysicalModel.Gravity, Zone);
                    uvws.AddValueProvider(value3D1Base);
                }
                else
                {
                    LxConsole.Error("Boussinesq模型要求密度为常数");
                    ComputeState.ReportInterrupt_FindError();
                }
            }
            var vis = Material.Viscosity.GetCopy(Zone);
            Value1Base<double> UVWgamma = GammaCommon == null ? vis : GammaCommon;
            Velocity.SetZone_CommonTerm(UVWgamma, Material.Viscosity, uvws.MakeValueList(), PhysicalModel.InitialVelocity.GetCopy(Zone), 0.5);
            Velocity.SetZone_ConvTerm(Density, PhysicalModel.ConvectionScheme);
            Velocity.SetZone_TransTerm(AppRes.SolveCenter.IsTransient, Value1Const.Double1, Density);
            if (AppRes.SolveCenter.IsTransient)
            {
                AppRes.SolveCenter.InitializeGroup.Initialization_Transient9.Add(Velocity.U.TransData);
                AppRes.SolveCenter.InitializeGroup.Initialization_Transient9.Add(Velocity.V.TransData);
                AppRes.SolveCenter.InitializeGroup.Initialization_Transient9.Add(Velocity.W.TransData);
                AppRes.SolveCenter.InitializeGroup.Initialization_Transient9.Add(Velocity.TransVelData);
            }

            PCor.SetZone(0.5, PhysicalModel.SolvingMethod, Velocity, P, Density);
            P.SetZone(new Value1Const.Double(0), Velocity, PCor, Density);
            Velocity.SetZone_RhieChow_AfterSetTrans(PCor, P);


            foreach (var t in Zone.BoundThreads)
            {
                switch (t.BoundaryType)
                {
                    case BoundaryTypes.FluidWall:
                        if (t.IsNoSlipingWall)
                        {
                            if (t.IsStaticWall)
                            {
                                Velocity.SetBoundary_NoSlipStaticWall(t);
                            }
                            else
                            {
                                Velocity.SetBoundary_NoSlipMovingWall(t.Velocity, t);
                            }
                        }
                        else
                        {
                            Velocity.SetBoundary_SlipWall(t);
                        }
                        //PCor.SetBound_ChangableP(t);
                        break;
                    case BoundaryTypes.Symmetry:
                        Velocity.SetBoundary_SlipWall(t);
                        //PCor.SetBound_ChangableP(t);
                        break;
                    case BoundaryTypes.VelocityInlet:
                        if (t.IsNormalToBound)
                        {
                            Velocity.SetBound_VelocityInlet(t.VelocityMagnitude, t);
                        }
                        else
                        {
                            Velocity.SetBound_VelocityInlet(t.Velocity, t);
                        }
                        //PCor.SetBound_ChangableP(t);
                        break;
                    case BoundaryTypes.PressureInOut:
                        if (t.IsNormalToBound)
                        {
                            Velocity.SetBound_PressureInOut(t.Pressure, t);
                        }
                        else
                        {
                            Velocity.SetBound_PressureInOut(t.Pressure, t.Velocity, t);
                        }
                        //PCor.SetBound_FixP(t);
                        break;
                    default:
                        LxConsole.Error("未知边界类型");
                        break;
                }
            }
            PCor.SetBound_AfterSettingVelocityVector();
            SetParameter_TurbBase();
        }

        public class BoussinesqSource : Value1Base<Vector3d>
        {
            public double RefT;
            public double[] ScalarT;
            public Vector3d Total;
            public BoussinesqSource(double RefT, double ThermExpanCoeff, double RefDensity, Vector3d Gravity, FvZone zone)
            {
                this.RefT = RefT;
                this.Total = RefDensity * ThermExpanCoeff * Gravity;
                this.ScalarT = zone.CFDFields.GetArray(ScalarTypes.Temperature);
                if (ScalarT == null)
                {
                    LxConsole.Error(string.Format($"区域{zone.name}设置Boussinesq模型出错，未获取温度场"));
                }
            }
            public override Vector3d GetCellValue(int cellIndex)
            {
                return (RefT - ScalarT[cellIndex]) * Total;
            }
        }
    }
}

//var threads = zone.BoundThreads;
//for (int i1 = 0; i1 < threads.Count; i1++)

//SolverControl.SIMPLEGroup.AddVelocity(zone, zone.AMGSolver, velocityVector);
//SolverControl.SIMPLEGroup.AddPcor(zone, zone.AMGSolver, PCor);
//SolverControl.SIMPLEGroup.AddP(zone, zone.AMGSolver, P);
//public override void initialize_ifNeeded(int taskIndex)
//{
//    velocityVector.initialize_ifNeeded(taskIndex);
//    P.initialize_ifNeeded(taskIndex);
//    T?.initialize_ifNeeded(taskIndex);
//}
//public override void initialize_force(int taskIndex)
//{
//    velocityVector.initialize_force(taskIndex);
//    P.initialize_force(taskIndex);
//    T?.initialize_force(taskIndex);
//}
//public override void updateMatrix(int taskIndex)
//{

//}
//public override void solve(int taskIndex)
//{
//    //1.求解动量方程
//    velocityVector.SetUVW_x(taskIndex);
//    velocityVector.U.CalInnerGradPhi(taskIndex);
//    velocityVector.V.CalInnerGradPhi(taskIndex);
//    velocityVector.W.CalInnerGradPhi(taskIndex);
//    velocityVector.SetUVW_b(taskIndex);
//    velocityVector.AddUVW_b_ViscositySource(taskIndex);
//    solveUVW(velocityVector.U, taskIndex);
//    solveUVW(velocityVector.V, taskIndex);
//    solveUVW(velocityVector.W, taskIndex);
//    //2.Rhie-Chow求解massImbalance
//    velocityVector.CalFaceMassFlux(taskIndex);
//    velocityVector.CalMassImbalance(taskIndex);
//    //3.求解压力修正
//    PCor.set_x_Sync(taskIndex);
//    PCor.set_b_Sync(taskIndex);
//    PCor.set_A_LU_Sync(taskIndex);
//    PCor.set_A_Diag_Sync(taskIndex);
//    PCor.cal_x_GS_Sync(taskIndex);
//    PCor.cal_x_GS_Sync(taskIndex);
//    PCor.cal_x_GS_Sync(taskIndex);
//    TasksControl.syncTasksWithOneAction(taskIndex, print);
//    //4.进行修正
//    PCor.predict_bound_Sync(taskIndex);
//    PCor.CalInnerGradPhi(taskIndex);
//    PCor.applyPressureCorrection_Sync(taskIndex);
//    //5.压力场,流场更新
//    P.Predict_bound(taskIndex);
//    P.CalInnerGradPhi(taskIndex);
//    velocityVector.UpdateUVWBound(taskIndex);
//}
//private void solveUVW(ScalarOfSteadyFlow UVW, int taskIndex)
//{
//    UVW.Set_A_LU_DiffPart(taskIndex);
//    UVW.set_A_Diag_Sync(taskIndex);
//    UVW.CorrectMassImbalance(taskIndex);
//    UVW.CalFaceNonOrthogonal(taskIndex);
//    UVW.CorrectNonOrthogonal(taskIndex);
//    UVW.Relax(taskIndex);
//    UVW.Cal_x_GS(taskIndex);
//}
//private void print()
//{
//    FindError.FindMaxAndMin(velocityVector.massImbalance);
//    FindError.FindMaxAndMin(PCor.x);
//}

//private void simpleAlg_test1_part2(int taskIndex)
//{
//    //1.求解动量方程
//    //SIMPLE.setUVW_x_All(taskIndex);
//    //U.calInnerGradPhi_Sync(taskIndex);
//    //V.calInnerGradPhi_Sync(taskIndex);
//    //W.calInnerGradPhi_Sync(taskIndex);
//    //SIMPLE.setUVW_b_All(taskIndex);
//    //SIMPLE.addUVW_b_ViscositySource_All(taskIndex);
//    //U.set_A_LU_DiffPart_Sync(taskIndex);
//    //U.set_A_Diag_Sync(taskIndex);
//    //U.correctMassImbalance_Sync(taskIndex);
//    //U.calFaceNonOrthogonal_Sync(taskIndex);
//    //U.correctNonOrthogonal_Sync(taskIndex);
//    //U.relax_Sync(taskIndex);
//    //// Console.WriteLine("U199计算值：{0}",(U.b[199] - ownerZone.cells[199].calALUProductValue(ownerZone.faces, U.A_LU_O,U. A_LU_N, U.x)) / U.A_Diag[199]);
//    //U.cal_x_GS_Sync(taskIndex);


//    //V.set_A_LU_DiffPart_Sync(taskIndex);
//    //V.set_A_Diag_Sync(taskIndex);
//    //V.correctMassImbalance_Sync(taskIndex);
//    //V.calFaceNonOrthogonal_Sync(taskIndex);
//    //V.correctNonOrthogonal_Sync(taskIndex);
//    //V.relax_Sync(taskIndex);
//    //V.cal_x_GS_Sync(taskIndex);

//    //W.set_A_LU_DiffPart_Sync(taskIndex);
//    //W.set_A_Diag_Sync(taskIndex);
//    //W.correctMassImbalance_Sync(taskIndex);
//    //W.calFaceNonOrthogonal_Sync(taskIndex);
//    //W.correctNonOrthogonal_Sync(taskIndex);
//    //W.relax_Sync(taskIndex);
//    //W.cal_x_GS_Sync(taskIndex);

//    //2.Rhie-Chow求解massImbalance
//    //SIMPLE.calFaceMassFlux_All(taskIndex);
//    //SIMPLE.calMassImbalance_All(taskIndex);
//    //Console.WriteLine("MassImbanlance :");
//    //FindError.FindMaxAndMin(SIMPLE.massImbalance);
//    //3.求解压力修正
//    //SIMPLE.set_x_All(taskIndex);
//    //SIMPLE.set_b_All(taskIndex);
//    //SIMPLE.set_A_LU_All(taskIndex);
//    //SIMPLE.set_A_Diag_All(taskIndex);
//    //SIMPLE.cal_x_GS_All(taskIndex);
//    //SIMPLE.cal_x_GS_All(taskIndex);
//    //SIMPLE.cal_x_GS_All(taskIndex);
//    //Console.WriteLine("PressureCorrection :");
//    //FindError.FindMaxAndMin(SIMPLE.x);

//    ////4.进行修正
//    //SIMPLE.calInnerGradPressureCorrection(taskIndex);
//    //SIMPLE.addPressureCorrection_All(taskIndex);
//    ////Console.WriteLine("Pressure :");
//    ////MaxAndMin.Write(SIMPLE.Pressure);
//    ////5.压力场,流场更新
//    //SIMPLE.updatePressureBound_All(taskIndex);
//    //SIMPLE.calGradPressure(taskIndex);
//    //SIMPLE.updateUVWBound_All(taskIndex);

//}
//private void simpleAlg_test2_part2()
//{
//    //1.求解动量方程
//    SIMPLE.setUVW_x_All();
//    U.calInnerGradField();
//    V.calInnerGradField();
//    W.calInnerGradField();
//    SIMPLE.setUVW_b_All();
//    // Console.WriteLine("1.U_b:max{0},min{1}", U.b.Max(), U.b.Min());
//    SIMPLE.addUVW_b_ViscositySource_All();
//    // Console.WriteLine("2.U_b:max{0},min{1}",U.b.Max(),U.b.Min());
//    //Console.WriteLine("b:{0}", U.b[210]);
//    U.set_A_LU_DiffPart_All();
//    //for(int i1 = 644; i1 <= 647; i1++)
//    //{
//    //    Console.WriteLine("alu o diff:{0},massflux:{1}", U.A_LU_O[i1], SIMPLE.massFlux[i1]);

//    //}

//    U.set_A_Diag_All();
//    //Console.WriteLine("adiag:{0}", U.A_Diag[210]);
//    //Console.WriteLine("b:{0}", U.b[210]);
//    //for (int i1 = 644; i1 <= 647; i1++)
//    //{
//    //    Console.WriteLine("alu o:{0}", U.A_LU_O[i1]);
//    //}
//    //Console.WriteLine("x=:{0}",(U.b[210]-U.A_LU_O[644]- U.A_LU_O[645]- U.A_LU_O[646]- U.A_LU_O[647])/U.A_Diag[210]);
//    V.set_A_LU_DiffPart_All();
//    V.set_A_Diag_All();
//    W.set_A_LU_DiffPart_All();
//    W.set_A_Diag_All();
//   // SIMPLE.addUVW_A_Diag_b_WallShearStress_All();
//   //foreach (int i in ownerZone.cells[210].CellFacesAsO)
//   // {
//   //     Console.WriteLine("x Index:{0},value:{1}", ownerZone.faces[i].NeighborCellIndex,U.x[ownerZone.faces[i].NeighborCellIndex]) ;
//   // }
//    U.correctMassImbalance_All();
//    //Console.WriteLine("质量不平衡修正后x=:{0}", (U.b[210] - U.A_LU_O[644] - U.A_LU_O[645] - U.A_LU_O[646] - U.A_LU_O[647]) / U.A_Diag[210]);
//    U.calFaceNonOrthogonal_All();
//    // Console.WriteLine("faceNonOrth:max:{0},min:{1}",ownerZone.faceCommonArray.Max(), ownerZone.faceCommonArray.Min());
//    U.correctNonOrthogonal_All();
//   // Console.WriteLine("非正交修正后x=:{0}", (U.b[210] - U.A_LU_O[644] - U.A_LU_O[645] - U.A_LU_O[646] - U.A_LU_O[647]) / U.A_Diag[210]);
//    U.relax_All();
//    U.cal_x_Jacobi_All_test();
//    //Console.WriteLine("计算出来的x210:{0}", U.x[210]);
//    //Console.WriteLine("计算出来的统计U :max:{0},min:{1},ave:{2}",U.x.Max(),U.x.Min(), U.x.Average());
//    //double maxU=0;
//    //int maxUCell=0;
//    //for(int i1 = 0; i1 < ownerZone.cells.Length; i1++)
//    //{
//    //    if (U.x[i1] > maxU)
//    //    {
//    //        maxUCell = i1;
//    //        maxU = U.x[i1];
//    //    }
//    //}
//    //Console.WriteLine("maxU:{0},位于：{1}",maxU,maxUCell);
//    //foreach (var t in ownerZone.BoundThreads)
//    //{
//    //    double bfa = 0;
//    //    for (int i1 = t.ThreadIndex.StartFaceIndex; i1 <= t.ThreadIndex.EndFaceIndex; i1++)
//    //    {
//    //        bfa += U.x[ownerZone.faces[i1].OwnerCellIndex];
//    //    }
//    //    Console.WriteLine(t.name + ":sum:{0}", bfa);
//    //}
//    V.correctMassImbalance_All();
//    V.calFaceNonOrthogonal_All();
//    V.correctNonOrthogonal_All();
//    V.relax_All();
//    V.cal_x_Jacobi_All_test();

//    W.correctMassImbalance_All();
//    W.calFaceNonOrthogonal_All();
//    W.correctNonOrthogonal_All();
//    W.relax_All();
//    W.cal_x_Jacobi_All_test();

//    //Console.WriteLine(U.x.Max());
//    //int errorcount = 0;
//    //foreach (var mx in U.x)
//    //{
//    //    if (mx > 2)
//    //    {
//    //        errorcount++;
//    //    }
//    //}
//    //Console.WriteLine(errorcount);
//    //2.Rhie-Chow求解massImbalance
//    SIMPLE.calMassFlux_temp();
//    SIMPLE.calMassImbalance_All();
//    ////3.求解压力修正
//    //SIMPLE.set_x_All();
//    //SIMPLE.set_b_All();
//    //SIMPLE.set_A_LU_All();
//    //SIMPLE.set_A_Diag_All();
//    //SIMPLE.cal_x_GS_All();
//    ////4.进行修正
//    //SIMPLE.addPressureCorrection_All();
//    ////5.压力场,流场更新
//    //SIMPLE.updatePressureBound_All();
//    //SIMPLE.calInnerGradPressure();
//    SIMPLE.updateUVWBound_All();
//    //Console.WriteLine(U.x.Max());
//    //errorcount = 0;
//    //foreach (var mx in U.x)
//    //{
//    //    if (mx > 2)
//    //    {
//    //        errorcount++;
//    //    }
//    //}
//    //Console.WriteLine(errorcount);
//}
//private void testFlowMatrixInitialize()
//{
//    SIMPLE.setZone_Density(new ValueConst(1));
//    SIMPLE.setZone_UVW(new ValueConst(10), new Value3dConst(new Vector3d()), new Value3dConst(new Vector3d(0, 0, 5)), 0.7, ownerZone.continuum.convectionScheme);
//    SIMPLE.initializeUVW_All();
//    SIMPLE.calMassFlux_temp();
//    SIMPLE.calMassImbalance_All();

//    T = new ScalarOfSteadyFlow();
//    T.createArray(ownerZone);
//    T.setMassDataArray(SIMPLE.massFlux, SIMPLE.massImbalance,SIMPLE.VelocityDirectionCorrection);
//    T.setZone(new ValueConst(10), new ValueConst(0), new ValueConst(300), 0.9, ownerZone.continuum.convectionScheme);
//    var threads = ownerZone.BoundThreads;
//    for (int i1 = 0; i1 < threads.Count; i1++)
//    {
//        switch (threads[i1].boundaryType)
//        {
//            case BoundaryTypes.流体壁面:
//                ValueConst fixT = new ValueConst(threads[i1].T);
//                T.setBoundary_fixValue(fixT, i1, threads[i1]);
//                break;
//            case BoundaryTypes.对称:
//                T.setBoundary_symmetry(i1, threads[i1]);
//                break;
//            default:
//                break;
//        }
//    }
//    T.initialize_All();
//}
//private void testFLowMatrixSolve()
//{
//    T.set_x_All();
//    T.calInnerGradField();
//    T.set_b_All();
//    T.set_A_LU_DiffPart_All();
//    T.set_A_Diag_All();
//    T.correctMassImbalance_All();
//    T.calFaceNonOrthogonal_All();
//    T.correctNonOrthogonal_All();
//    T.relax_All();
//    ownerZone.MultiGridSolver.Solve_VCycle(T);
//    T.cal_x_GS_All();
//    T.cal_x_GS_All();
//    T.updateBound_All();
//}